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Summary 

A viscoelastic higher-order thick beam finite element formulation is extended to 
include elastodynamic deformations. The material constitutive law is a special 
differential form of the Maxwell solid. In the constitutive model, the elastic 
strains and the conjugate viscous strains are coupled through a system of first- 
order ordinary differential equations. The total time-dependent stress is the 
superposition of its elastic and viscous components. The elastodynamic 
equations of motion are derived from the virtual work principle. Computational 
examples are carried out for a thick orthotropic cantilevered beam. A quasi-static 
relaxation problem is employed as a validation test for the elastodynamic 
algorithm. The elastodynamic code is demonstrated by analyzing the damped 
vibrations of the beam which is deformed and then released to freely vibrate. 

Introduction 

The combination of highly viscous low-modulus materials with the traditional 
high-modulus materials produces stiff, highly damped load carrying structures. 
The quasi-static and dynamic analyses of such structures require improvements in 
the material damping representation over the velocity proportional damping 
schemes. Halpin and Pagano 1 demonstrated that the relaxation moduli for 
anisotropic solids produce symmetric matrices that can be expanded in a Prony 
series form (i.e., a series of exponentially decaying terms). Early viscoelastic 
models for small deformations of composites focused on computing the complex 
moduli for anisotropic solids from the elastic properties of the fibers and the 
complex modulus properties of the matrix. 2,3 Recently, various classical 
constitutive models have been used including generalized Maxwell and Kelvin- 
Voigt solids. 4,5 These constitutive models have practical value since they provide 
adequate approximations for the dynamic softening and hysteresis effects - the 
phenomena that are not directly proportional to strain rates. 

Coleman and Noll 6 and Schapeiy 7 presented comprehensive discussions on the 
history integral form of the Maxwell solid. Johnson and Stacer 8 developed a 
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differential form of the Maxwell solid constitutive law for large strain viscoelastic 
deformations of rubber. The formulation required additional viscous strain 
variables that are conjugate to the elastic strains. The same constitutive law was 
used by Johnson et al. 9,10 to formulate a viscoelastic, large-displacement shell 
finite element. In the work of Johnson and Tessler, 11 the differential constitutive 
law was implemented within Tessler’s 12 higher-order beam theory, giving rise to a 
quasi-static viscoelastic thick beam finite element formulation. They presented 
several quasi-static numerical solutions demonstrating the relaxation, creep, and 
cyclic creep response of thick beams. They demonstrated that only minor 
modifications to an elastic finite element code are needed to produce a 
computationally efficient viscoelastic code. 

In this paper, the quasi-static viscoelastic formulation of Johnson and Tessler 11 
is extended to elastodynamics. The mass matrix for the higher-order beam is 
derived and a modified Newmark algorithm is employed to integrate the equations 
of motion. The modified Newmark algorithm contains an implicit trapezoidal 
integration scheme which enables accurate integration of the first-order differential 
equations for the viscous strains. 

Viscoelastic Higher-Order Beam 

In this section, an internal variable formulation for a Maxwell solid higher-order 
beam theory 11 is outlined. The beam dimensions and sign convention are shown in 
Figure 1 . The differential form of the viscoelastic constitutive model for the beam 
is written in matrix form as 


( 1 ) 
( 2 ) 
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The vectors s and e are the engineering stresses and strains. The vectors *e 
(n = 1,2,..., N) are the internal variables, i.e., conjugate viscous strains, where N 
is the number of terms in a Prony series representation of the material’s stress 
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relaxation response. The matrices C and *C contain the elastic and viscous 
stiffness coefficients. In the higher-order beam theory, 12 the components of the 
displacement vector are approximated through the beam thickness by way of five 
kinematic variables, i.e., 


u x (x,z,t) = u(x,t) + h£d(x,t) 


( 


•2 I'l 


U z (x,z,t) = w(x,t) + £w x (x,t)+\£‘ --jw 2 (x,t) 


( 3 ) 

( 4 ) 


where £ = z/h denotes a nondimensional thickness coordinate and 2h is the total 
thickness. The function u{x,t) represents the midplane (i.e. reference plane) axial 
displacement, 6{x,t) is the bending rotation of the cross-section of the beam, 
w(x,t) is the weighted-average deflection, and and w 2 (jc, t) are the higher- 

order transverse displacement variables enabling a parabolic distribution of 
u z (x,z,t ) through the thickness. 



Figure 1. Thick-beam geometry, kinematics, and loading. 


The above displacement assumptions give rise to the following thickness 
distributions for the strains: a linear axial strain, a cubic transverse normal strain, 
and a quadratic transverse shear strain. These strain components have the 
following form 
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where ^(0 = hv B f(4-7f 2 )/l7, 0,(0 = 140(3 -f 2 )/l7, fe<D = 5(l-f 2 )/4, 

and v 13 is Poisson’s ratio. The simplest finite element approximation of this 
beam theory involves a three-node configuration (see Figure 2) which is achieved 
by the following interpolations 

u{r],t) = (1 - tj )u e 0 (t) + 77 iif ( 0 , 6 /( 17 , t) = (1 - 77 ^( 7 ) + 776 // (7), 

w(T 7 , 7) = (1 - 77 ) Wq (7) + 77w/ ( 7 ) - 1 77(1 - r7)(^o (0 ~ (0) . ( 8 ) 

w 1 (t 7 ,7) = W/(7), w 2 (77,7) = W{(t) 

where r\=-x!l is the nondimensional axial coordinate. The nodal degrees-of- 
freedom at the two ends of the element are subscripted with indices 0 and 1. Since 
the strains do not possess derivatives of the Wj ( 77 , 7 ) and vv 2 ( 77 , 7 ) variables, these 
variables need not be continuous at the element nodes and, hence, their simplest 
approximation is constant for each element. Their corresponding degrees-of- 
freedom are attributed to a node at the element midspan. 



For dynamic loading, the virtual work statement for an element of volume V 
with the differential Maxwell constitutive law can be written as 

N 

| p(u x 8u x + u z 8u z ) dV + 1 e T CSe dV + if: e T *CS* n edV-8W = 0 (9) 

n - 1 
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where the first integral represents the virtual work done by inertial forces, the 
second is the internal virtual work done by the elastic stresses, the third is the 
internal virtual work done by the viscous stresses, and 8W is the virtual work 
done by the external forces. Introducing ( 8 ) into (3)-(4) and substituting the 
results into (5)-(7) yields finite element approximations of the strains in terms of 
the nodal variables, i.e., 
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( 10 ) 


( 11 ) 


and u r = (u o ,w o ,0 o ,W l ,W 2 ,u l ,w i ,dy) e denotes the element nodal displacement 
vector. Next, a set of analogous nodal variables, *u, and corresponding viscous 
strains, *e, are introduced. These are related by 


* 

n 


e = B!u 


( 12 ) 


The *u variables, which carry the time dependent information for the material 
within the element, are independent from element to element. The displacements 
u x and u z are then expressed in terms of the element nodal degrees-of-freedom 
using (3), (4) and ( 8 ), giving rise to u x = <£ x t ii and u z = O z T u, where < 3 > x (£, 77 ) and 
O z (^, 77 ) are vectors of the interpolation functions. The virtual work statement 
for an element then becomes 

ii T J p(O x <D x r + O.O, 7 ) dV da + u r J B r C B dV Sa 

n , (13) 

+ ^* n u T jB T *CBdV8* n u-8W=0 

n = 1 


By defining the integrals in (13) as the mass, m, elastic stiffness, k, and viscous 
stiffness, *k, matrices, there results 
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|u r m + u T k]&i + 
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5> r *k 


.«= 1 


S:u-8W = 0 


(14) 


Since <5u =<5*u when t is constant, the virtual work takes on a simpler form 


N 


ii r m + u 7 k + ]£ *u r *k 

n - 1 


3a-SW = 0 


(15) 


This implies that at any time t the element equilibrium equations are 

N 


mii + ku = f-^’k *u 

n = 1 


for each element 


(16) 


where f denotes the element consistent load vector due to the external loading. 
Introducing (10) and (12) into the differential equations for the strain variables in 
(2) yields 

*u + *u/r„ = u for each n (17) 

The global equilibrium equations are determined by the standard assembly of 
the element equations, (16). Note, there is no assembly for (17). The global 
equations of motion are 

Mii g + Ku g = F mech - F visc (18) 

where denotes the global nodal variable vector, M is the mass matrix, K is 
the elastic stiffness matrix, F mech is the global force vector due to mechanical 

N 

loads, and F Ww . is the assembled vector for *k *u. The viscoelastic problem is 

n = 1 

solved by simultaneously integrating the first order differential equations, (17), 
and the second order equations, (18), where the latter is subject to the appropriate 
boundary restraints. 

As far as the finite element implementation is concerned, a conventional linear 
elastic code can be readily adapted to perform a dynamic analysis for a structure 
made from a Maxwell material, i.e., a material whose relaxation stiffness 
coefficients can be modeled with a Prony series. First, the viscous stiffness 
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coefficients, *C lj , are used in place of the elastic values to compute the element 

viscous stiffness matrices, *k, which are stored for repeated use. The internal 
nodal variables for each element, *u, are set equal to their initial values. A 
modified Newmark algorithm is then used to integrate (18). The modification is 
required so that (17), stiff relaxation equations, are implicitly integrated. 

Applications 

Two numerical solutions representative of quasi-static and free vibration 
deformations of a cantilevered thick orthotropic beam are presented. The beam 
dimensions are as follows: L = 0.1m, 2 h = 0.02m, and b = 1.0m (refer to Figure 
1). The beam elastic stiffness coefficients for the state of plane stress can be 
written in terms of engineering material constants as 

Qi = Ej(\ - v^v^), C 33 = E z /( 1 - v^Vy.), C J3 = v Jtz C 33 , C 55 = G xz 

A unidirectional E-glass/epoxy laminate is considered for which the material 
constants are: E x = 38.6 GPa, E z = 8.27 GPa, G xz = 4.14 GPa, v xz = 0.26, and 
p = 1.8 g/cm 3 . The viscous relaxation properties are computed from complex 
modulus versus frequency data for the E-glass/epoxy. 13 The equations for the real 
and imaginary parts of the effective modulus for a ten term Maxwell solid 14 are 
least-squares fit to the measured data in a frequency range of 45 Hz - 145 Hz. 
The least squares minimization was performed using a quadratic programming 
method that enforced the constraint that each of the moduli in the Maxwell solid 
must be positive. The time constants for the ten terms are; t„= 10 4 , 10' 3 5 , 10' 3 , 
10' 2 ' 5 , 10' 2 , ... , 10 1 , and infinity. The resulting Maxwell solid had its moduli 
scaled so that its equivalent nondimensional Prony series has the form 

P(t) = 1.0 + 0.01755 e" 0 0001 ' +0.000257 e~°* u +0.072014 e ~°- 3162 ‘ 

The time dependent stiffness values for the E-glass/epoxy are then given by 
*C = P(t)C. 

Example 1 . A cantilever beam (reference Figure 1) with w, u, 0 fixed at point A 
has a prescribed deflection W at point B that is ramped from 0 to -1 cm in 0.05 
sec and then held constant at - 1 cm. Figure 3 depicts the maximum axial stress at 
point C as a function of time. Also shown are the elastic and viscous stress 
components comprising the total stress. The decay of the total viscoelastic stress 
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to its elastic value as time is increased demonstrates the expected step-strain 
relaxation behavior. Note, when this problem is solved with a velocity 
proportional viscous damping model the viscous stress for time greater than 0.05 
sec is zero. 



Time (sec) 

Figure 3. Cantilevered beam under prescribed tip deflection. 

Stress in top surface at support (point C). 

Example 2 . The cantilevered beam of Example 1 has its tip (point B) released at t 
= 0.1 sec and is allowed to vibrate freely thereafter. The resulting tip deflection 
as a function of time is shown in Figure 4. Figure 5 depicts the value of the 
dynamic axial stress at point C as a function of time. Note that both relaxation 
and damped vibration response are obtained within the same finite element 
solution. 


Conclusions 

An elastodynamic formulation, which includes a differential form of the 
Maxwell viscous solid constitutive theory, has been implemented in a higher- 
order-theory beam finite element. The attractive features of the formulation 
include: (1) The constitutive constants are the same as those needed in the 
classical history-integral model, and they are also readily available from step-strain 
relaxation tests, (2) The internal variables are conjugate to the elastic strain 
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measures; hence, they are consistent with the kinematic assumptions of the elastic 
formulation, (3) The update of the state variables can be performed in a parallel 
computing environment, allowing the viscous force vector in the equations of 
motion to be determined efficiently within the modified Newmark algorithm, 


Number of elements = 40 Time Step - 1.0e-5 seconds 
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Figure 4. Cantilevered beam under prescribed tip deflection. 
Tip released - vibration of tip (Point B). 
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Figure 5. Cantilevered beam under prescribed tip deflection. 

Tip released - Stress on top surface at support. (Point C). 
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(4) Applications of time-dependent displacements and loads are performed within 
the same finite element algorithm, and (5) The higher-order beam theory accounts 
for both transverse shear and transverse normal deformations — the effects that 
need to be accounted for in thick and highly orthotropic beams, and in high- 
frequency dynamics. 

Example 2, in which a quasi-static enforced deformation is followed by a high 
frequency free vibration, demonstrates the capability of this finite element 
formulation to simulate physically important phenomena that is computationally 
difficult to obtain. 
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